      program ramix2d
      implicit real*8 (a-h,o-z)
      real*8 lb
      logical max
      parameter(np=4,mp=np+1,nplay=101,ndec=12,
     &          ftol=1.0d-12,neps=4)
c
      external famoeb
      character*6 name(np)
      character*24 fdate
      dimension p(mp,np),x(np),y(mp)
      dimension is(nplay,ndec),pce(nplay)
      dimension lb(np),ub(np),c(np),vm(np),xopt(np),xp(np),
     &   fstar(neps),nacp(np)
      common /mats/is
      common /pchm/pce
c
      data name/'BETA','GAM','A0','A1'/
c
      open(unit=10,file='ra2011.data')   !with card sleeves on boards
c
      read(10,*)((is(i,n),n=1,ndec),i=1,nplay)
c
c      do 5 id=1,nplay
c         write(*,'(i2,2x,12i3)')id,(is(id,j),j=1,ndec)
c    5 continue
c      stop
c
c---------------------------------------------------------
      write(*,'('' ***** ramix2d *****'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/"   All 2011 Sessions")')
c
      call fpeu(pce)
c
c Set parameters for simulated annealing
      max = .true.
      eps = 1.d0
      rt = 0.5
      iseed1 = 8069
      iseed2 = 4012
      ns = 20
      nt = 25
      maxevl = (nt+5)*nt*ns*np/10 + 2
c      maxevl = 1
      iprint = 1
      do 30 i = 1, np
         lb(i) = -5.d0
         ub(i) =  5.d0
         c(i) = 2.0
   30 continue
c initial values for sa
         x(1) = dlog(0.9d0/0.7d0 - 1.d0)
         x(2) = dlog(2.d0/1.d0 - 1.d0) 
         x(3) = dlog(1.d0/0.5d0 - 1.d0)
         x(4) = dlog(0.5d0/0.1d0 - 1.d0)
c
      T = 8.d0
      do 40, i = 1, np
         vm(i) = 0.5d0*(ub(i)-lb(i))
   40 continue
      write(*,1000) t, rt, eps, ns, nt, neps, maxevl
c
      call sa(np,x,max,rt,eps,ns,nt,neps,maxevl,lb,ub,c,iprint,iseed1,
     &        iseed2,t,vm,xopt,fopt,nacc,nfcnev,nobds,ier,
     &        fstar,xp,nacp)
c
c      write(*,'(/,''  ****   results after sa   ****   '')') 
      call prtvec(xopt,np,'solution')
      call prtvec(vm,np,'final step length')
      write(*,'(/" optimal function value: ",g15.7)')fopt
c
1000  format(/,' simulated annealing summary',/,
     &       /,' initial temp: ', g8.2, '   rt: ',g8.2, '   eps: ',g8.2,
     &       /,' ns: ',i3, '   nt: ',i2, '   neps: ',i2,' maxevl: ',i10)
c
c reverse parameter transformation
         do 45 k=1,np
            x(k)=xopt(k)
   45    continue
c
c initialize starting vertices of the simplex. 
         call init(x,p,y,famoeb)
c 
c call the simplex routine. 
         ndim=np 
         iter=3000
c         iter=1
         call amoeba(p,y,mp,np,ndim,ftol,famoeb,iter) 
c 
c output results and stop. 
         write(*,'(a,i5)')'iterations: ',iter 
c compute average values of vertices. 
         call avg(p,x,np,mp) 
         xllf=-famoeb(x) 
c reverse parameter transformation
         x(1)=0.9d0/(1.d0+dexp(x(1)))
         x(2)=2.8d0/(1.d0+dexp(x(2)))+2.2d0
         x(3)=1.d0/(1.d0+dexp(x(3)))
         x(4)=(1.d0-x(3))/(1.d0+dexp(x(4)))
         a2=1.d0-x(3)-x(4)
c
      write(*,'(//"These are the sa-amoeba parameter estimates:"/)')
      do 50 k=1,np
         write(*,'(2x,a6,2x,g13.5)')name(k),x(k)
   50 continue
      write(*,'(2x,a6,2x,g13.5)')"A2    ",a2
      write(*,'(/"Maximized LL = ",g13.6)')xllf
c
      call xpost(x)
      end
c***********************************************************************
      double precision function famoeb(x)
      implicit real*8 (a-h,o-z)
      parameter (np=4)
      dimension x(np)
c
      beta=0.9d0/(1.d0+dexp(x(1)))
      gam=2.8d0/(1.d0+dexp(x(2)))+ 2.2d0
      a0=1.d0/(1.d0+dexp(x(3)))
      a1=(1.d0-a0)/(1.d0+dexp(x(4)))
c
      call fam(beta,gam,a0,a1,pt)

      famoeb = -pt
c
      return
      end
c***********************************************************************
      subroutine fam(beta,gam,a0,a1,pt)
      implicit real*8 (a-h,l,o-z)
      parameter (nplay=101,ndec=12)
      dimension pce(nplay),paa(nplay)
      common /pchm/pce

      call fpaa(beta,gam,paa)
      p0 = 0.5d0**ndec
      a2 = 1.d0-a0-a1
c
      pt = 0.d0
      do 50 id=1,nplay
         pid = a0*p0 + a1*paa(id) + a2*pce(id)
         pid = dlog(pid)
         pt = pt + pid
   50 continue
      return
      end
c***********************************************************************
      subroutine fpeu(pce)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,nplay=101)
      dimension z(2,ndec),pm(2,ndec),pce(nplay)
      dimension is(nplay,ndec)
      common /mats/is
c
      beta = 1.d0
      gam = 10.d0
      call ystar(beta,z)
      call fpm(gam,z,pm)
c      write(*,'(12g10.3)')(pm(1,i),i=1,ndec)
c      stop
c
      do 70 id=1,nplay
         pin = 1.d0
         do 60 j=1,ndec
            ix = is(id,j)
            pin = pin*pm(ix,j)
   60    continue
         pce(id) = pin
c      write(*,'(i2," pce = ",g13.5)')id,pce(id)
   70 continue
c      stop
      return
      end
c***********************************************************************
      subroutine fpaa(beta,gam,paa)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,nplay=101)
      dimension z(2,ndec),pm(2,ndec),paa(nplay)
      dimension is(nplay,ndec)
      common /mats/is
c
      call ystar(beta,z)
      call fpm(gam,z,pm)
c      write(*,'(12g10.3)')(pm(1,i),i=1,ndec)
c      stop
c
      do 70 id=1,nplay
         pin = 1.d0
         do 60 j=1,ndec
            ix = is(id,j)
            pin = pin*pm(ix,j)
   60    continue
         paa(id) = pin
c      write(*,'(i2," paa = ",g13.5)')id,paa(id)
   70 continue
c      stop
      return
      end
c***********************************************************************
      subroutine ystar(beta,z)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension z(2,ndec)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      z(1,1) = 0.5d0*y1
      z(2,1) = 0.5d0*beta*y1
      z(1,2) = 0.5d0*y1
      z(2,2) = 0.5d0*beta*y1
      z(1,3) = 0.5d0*y1
      z(2,3) = 0.5d0*beta*y2
      z(1,4) = 0.5d0*y1
      z(2,4) = 0.5d0*beta*y2
      z(1,5) = 0.5d0*y1
      z(2,5) = 0.5d0*beta*y3
      z(1,6) = 0.5d0*y1
      z(2,6) = 0.5d0*beta*y3
      z(1,7) = y1*beta*1.d0/3.d0
      z(2,7) = y1*1.d0/3.d0
      z(1,8) = y2*beta*1.d0/3.d0
      z(2,8) = y1*1.d0/3.d0
      z(1,9) = y3*beta*1.d0/3.d0
      z(2,9) = y1*1.d0/3.d0
      z(1,10) = y4*2.d0/3.d0
      z(2,10) = y4*pjb
      z(1,11) = y4*2.d0/3.d0
      z(2,11) = y5*pjb
      z(1,12) = y4*2.d0/3.d0
      z(2,12) = y6*pjb
c
      return
      end
c***********************************************************************
      subroutine fpm(gam,z,pm)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension z(2,ndec),y(2),pm(2,ndec)
c
      do 20 j=1,ndec
         y(1) = dexp(gam*z(1,j))
         y(2) = dexp(gam*z(2,j))
         den = y(1) + y(2)
         pm(1,j) = y(1)/den
         pm(2,j) = y(2)/den
   20 continue
c
      return
      end
c***********************************************************************
      subroutine xpost(x)
      implicit real*8 (a-h,l,o-z)
      parameter (np=4,nplay=101,ndec=12)
      dimension x(np),pce(nplay),paa(nplay)
      dimension pid(3),alf(3)
      dimension pe(2,ndec),pa(2,ndec),z(2,ndec),pst(ndec)
      common /pchm/pce
c
      beta=x(1)
      gam=x(2)
      a0=x(3)
      a1=x(4)
      a2=1.d0-a0-a1
      p0 = 0.5d0**ndec
      call fpaa(beta,gam,paa)
c
      call ystar(beta,z)
      call fpm(gam,z,pa)
      write(*,'(//"paa: ",12g12.5)')(pa(1,k),k=1,ndec)
      beta = 1.d0
      gam = 1.d1
      call ystar(beta,z)
      call fpm(gam,z,pe)
      write(*,'(/"peu: ",12g12.5)')(pe(1,k),k=1,ndec)
      do 5 j=1,ndec
         pst(j) = a0*0.5d0 + a1*pa(1,j) + a2*pe(1,j)
    5 continue
      write(*,'(/"pst: ",12g12.5)')(pst(k),k=1,ndec)
c
      write(*,'(//"Ex Post Individual:"/)')
      pt = 0.d0
      do 50 id=1,nplay
         pid(1) = a0*p0
         pid(2) = a1*paa(id)
         pid(3) = a2*pce(id)
         pit = pid(1) + pid(2) + pid(3)
         amax = -1.d0
         do 45 n=1,3
            alf(n) = pid(n)/pit
            if (alf(n).gt.amax) then 
               amax=alf(n)
               m=n-1
            endif
   45    continue
         pit = dlog(pit)
         pt = pt + pit
         write(*,'(i2,2x,4g13.4,2x,i3))')id,pit,(alf(n),n=1,3),m
   50 continue
      write(*,'(/"Total LL = ",g13.5)')pt
      return
      end
c***********************************************************************
      subroutine init(xint,p,y,famoeb)
      implicit real*8 (a-h,o-z) 
      parameter (np=4,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external famoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      do 2 k=1,np
c         if (xx(k).le.1.d-10) xx(k)=1.d-10
c         if (xx(k).ge.0.999) xx(k)=0.999
c         xxx(k)=dlog(1.d0/xx(k) - 1.d0)
c    2 continue
c
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
      y0=-famoeb(xxx)
      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
      delta=0.1d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=famoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c*********************************************************************
c ************** newlibrary of subroutines ***************************
      subroutine avg(p,x,np,mp)
      implicit real*8 (a-h,o-z)
      dimension p(mp,np),x(np)
c
c compute average values of vertices of the simplex.
      do 15 k=1,np
         do 14 j=2,mp
   14       p(1,k)=p(1,k)+p(j,k)
   15    x(k)=p(1,k)/mp
      return
      end
c------------------------------------------------------------------------
      subroutine amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*8 (a-h,o-z)
      parameter (nmax=20,alpha=1.d0,beta=0.5d0,gamma=2.d0)
      dimension p(mp,np),y(mp),pr(nmax),prr(nmax),pbar(nmax)
c
c This is the simplex algorithm from Numerical Recipes, with 
c modifications.
c When amoeba is called, iter must be set equal to the maximum
c number of iterations that will be allowed.  On return,
c iter equals the actual number of iterations.
c Thus, an example of a calling sequence would be:
c          itmax=5000
c          iter=5000
c          call aboeba(p,y,mp,np,ndim,ftol,funk,iter)
c          if (iter.ge.itmax) ...
c On return, if the if condition evaluates to true, convergence
c has not been acheived in the alloted number of iterations.
c     --PWW   1/26/94
c
c
      itmax=iter
      if (itmax.le.0) write(*,100)itmax
  100 format('*****WARNING: Maximum no. of iterations passed to',
     &       ' amoeba=',i6,'*****')
      mpts=ndim+1
      iter=0
    1 ilo=1
      if(y(1).gt.y(2))then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 11 i=1,mpts
        if(y(i).lt.y(ilo)) ilo=i
        if(y(i).gt.y(ihi))then
        inhi=ihi
        ihi=i
        else if(y(i).gt.y(inhi))then
        if(i.ne.ihi) inhi=i
        endif
   11 continue
      rtol=2.d0*dabs(y(ihi)-y(ilo))/(dabs(y(ihi))+dabs(y(ilo)))
      if(rtol.lt.ftol)return
c
c check number of iterations.
      if (iter.ge.itmax) then
         write(*,101)itmax
         return
         end if
  101 format('*****Maximum iterations (',i6,') exceeded in amoeba.'/
     &       '     Control passed to calling routine.*****')
      iter=iter+1
      do 12 j=1,ndim
        pbar(j)=0.d0
   12 continue
      do 14 i=1,mpts
        if(i.ne.ihi)then
        do 13 j=1,ndim
        pbar(j)=pbar(j)+p(i,j)
   13   continue
        endif
   14 continue
      do 15 j=1,ndim
        pbar(j)=pbar(j)/ndim
        pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
   15 continue
      ypr=funk(pr)
      if(ypr.le.y(ilo))then
        do 16 j=1,ndim
        prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
   16   continue
        yprr=funk(prr)
        if(yprr.lt.y(ilo))then
        do 17 j=1,ndim
        p(ihi,j)=prr(j)
   17   continue
        y(ihi)=yprr
        else
        do 18 j=1,ndim
        p(ihi,j)=pr(j)
   18   continue
        y(ihi)=ypr
        endif
      else if(ypr.ge.y(inhi))then
        if(ypr.lt.y(ihi))then
        do 19 j=1,ndim
        p(ihi,j)=pr(j)
   19   continue
        y(ihi)=ypr
        endif
        do 21 j=1,ndim
        prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
   21   continue
        yprr=funk(prr)
        if(yprr.lt.y(ihi))then
        do 22 j=1,ndim
        p(ihi,j)=prr(j)
   22   continue
        y(ihi)=yprr
        else
        do 24 i=1,mpts
        if(i.ne.ilo)then
        do 23 j=1,ndim
                pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=pr(j)
   23   continue
        y(i)=funk(pr)
        endif
   24   continue
        endif
      else
        do 25 j=1,ndim
        p(ihi,j)=pr(j)
   25   continue
        y(ihi)=ypr
      endif
      go to 1
      end
c******************************************************
      subroutine fcn(n,x,f)
      implicit real*8 (a-h,o-z)
      dimension x(n)
c
c recover the parameters
c
      beta=0.9d0/(1.d0+dexp(x(1)))
      gam=2.8d0/(1.d0+dexp(x(2))) + 2.2d0
      a0=1.d0/(1.d0+dexp(x(3)))
      a1=(1.d0-a0)/(1.d0+dexp(x(4)))
c
      call fam(beta,gam,a0,a1,f)
c
      return
      end
c***************************************************************************
C  Simulated Annealing Version: 3.2
C  Date: 1/22/94.
C  Differences compared to Version 2.0:
C     1. If a trial is out of bounds, a point is randomly selected
C        from LB(i) to UB(i). Unlike in version 2.0, this trial is
C        evaluated and is counted in acceptances and rejections.
C        All corresponding documentation was changed as well.
C  Differences compared to Version 3.0:
C     1. If VM(i) > (UB(i) - LB(i)), VM is set to UB(i) - LB(i).
C        The idea is that if T is high relative to LB & UB, most
C        points will be accepted, causing VM to rise. But, in this
C        situation, VM has little meaning; particularly if VM is
C        larger than the acceptable region. Setting VM to this size
C        still allows all parts of the allowable region to be selected.
C  Differences compared to Version 3.1:
C     1. Test made to see if the initial temperature is positive.
C     2. WRITE statements prettied up.
C     3. References to paper updated.
C
C  Synopsis:
C  This routine implements the continuous simulated annealing global
C  optimization algorithm described in Corana et al.'s article
C  "Minimizing Multimodal Functions of Continuous Variables with the
C  "Simulated Annealing" Algorithm" in the September 1987 (vol. 13,
C  no. 3, pp. 262-280) issue of the ACM Transactions on Mathematical
C  Software.
C
C  A very quick (perhaps too quick) overview of SA:
C     SA tries to find the global optimum of an N dimensional function.
C  It moves both up and downhill and as the optimization process
C  proceeds, it focuses on the most promising area.
C     To start, it randomly chooses a trial point within the step length
C  VM (a vector of length N) of the user selected starting point. The
C  function is evaluated at this trial point and its value is compared
C  to its value at the initial point.
C     In a maximization problem, all uphill moves are accepted and the
C  algorithm continues from that trial point. Downhill moves may be
C  accepted; the decision is made by the Metropolis criteria. It uses T
C  (temperature) and the size of the downhill move in a probabilistic
C  manner. The smaller T and the size of the downhill move are, the more
C  likely that move will be accepted. If the trial is accepted, the
C  algorithm moves on from that point. If it is rejected, another point
C  is chosen instead for a trial evaluation.
C     Each element of VM periodically adjusted so that half of all
C  function evaluations in that direction are accepted.
C     A fall in T is imposed upon the system with the RT variable by
C  T(i+1) = RT*T(i) where i is the ith iteration. Thus, as T declines,
C  downhill moves are less likely to be accepted and the percentage of
C  rejections rise. Given the scheme for the selection for VM, VM falls.
C  Thus, as T declines, VM falls and SA focuses upon the most promising
C  area for optimization.
C
C  The importance of the parameter T:
C     The parameter T is crucial in using SA successfully. It influences
C  VM, the step length over which the algorithm searches for optima. For
C  a small initial T, the step length may be too small; thus not enough
C  of the function might be evaluated to find the global optima. The user
C  should carefully examine VM in the intermediate output (set IPRINT =
C  1) to make sure that VM is appropriate. The relationship between the
C  initial temperature and the resulting step length is function
C  dependent.
C     To determine the starting temperature that is consistent with
C  optimizing a function, it is worthwhile to run a trial run first. Set
C  RT = 1.5 and T = 1.0. With RT > 1.0, the temperature increases and VM
C  rises as well. Then select the T that produces a large enough VM.
C
C  For modifications to the algorithm and many details on its use,
C  (particularly for econometric applications) see Goffe, Ferrier
C  and Rogers, "Global Optimization of Statistical Functions with
C  Simulated Annealing," Journal of Econometrics, vol. 60, no. 1/2, 
C  Jan./Feb. 1994, pp. 65-100.
C  For more information, contact 
C              Bill Goffe
C              Department of Economics and International Business
C              University of Southern Mississippi 
C              Hattiesburg, MS  39506-5072 
C              (601) 266-4484 (office)
C              (601) 266-4920 (fax)
C              bgoffe@whale.st.usm.edu (Internet)
C
C  As far as possible, the parameters here have the same name as in
C  the description of the algorithm on pp. 266-8 of Corana et al.
C
C  In this description, SP is single precision, DP is double precision,
C  INT is integer, L is logical and (N) denotes an array of length n.
C  Thus, DP(N) denotes a double precision array of length n.
C
C  Input Parameters:
C    Note: The suggested values generally come from Corana et al. To
C          drastically reduce runtime, see Goffe et al., pp. 90-1 for
C          suggestions on choosing the appropriate RT and NT.
C    N - Number of variables in the function to be optimized. (INT)
C    X - The starting values for the variables of the function to be
C        optimized. (DP(N))
C    MAX - Denotes whether the function should be maximized or
C          minimized. A true value denotes maximization while a false
C          value denotes minimization. Intermediate output (see IPRINT)
C          takes this into account. (L)
C    RT - The temperature reduction factor. The value suggested by
C         Corana et al. is .85. See Goffe et al. for more advice. (DP)
C    EPS - Error tolerance for termination. If the final function
C          values from the last neps temperatures differ from the
C          corresponding value at the current temperature by less than
C          EPS and the final function value at the current temperature
C          differs from the current optimal function value by less than
C          EPS, execution terminates and IER = 0 is returned. (EP)
C    NS - Number of cycles. After NS*N function evaluations, each
C         element of VM is adjusted so that approximately half of
C         all function evaluations are accepted. The suggested value
C         is 20. (INT)
C    NT - Number of iterations before temperature reduction. After
C         NT*NS*N function evaluations, temperature (T) is changed
C         by the factor RT. Value suggested by Corana et al. is
C         MAX(100, 5*N). See Goffe et al. for further advice. (INT)
C    NEPS - Number of final function values used to decide upon termi-
C           nation. See EPS. Suggested value is 4. (INT)
C    MAXEVL - The maximum number of function evaluations. If it is
C             exceeded, IER = 1. (INT)
C    LB - The lower bound for the allowable solution variables. (DP(N))
C    UB - The upper bound for the allowable solution variables. (DP(N))
C         If the algorithm chooses X(I) .LT. LB(I) or X(I) .GT. UB(I),
C         I = 1, N, a point is from inside is randomly selected. This
C         This focuses the algorithm on the region inside UB and LB.
C         Unless the user wishes to concentrate the search to a par-
C         ticular region, UB and LB should be set to very large positive
C         and negative values, respectively. Note that the starting
C         vector X should be inside this region. Also note that LB and
C         UB are fixed in position, while VM is centered on the last
C         accepted trial set of variables that optimizes the function.
C    C - Vector that controls the step length adjustment. The suggested
C        value for all elements is 2.0. (DP(N))
C    IPRINT - controls printing inside SA. (INT)
C             Values: 0 - Nothing printed.
C                     1 - Function value for the starting value and
C                         summary results before each temperature
C                         reduction. This includes the optimal
C                         function value found so far, the total
C                         number of moves (broken up into uphill,
C                         downhill, accepted and rejected), the
C                         number of out of bounds trials, the
C                         number of new optima found at this
C                         temperature, the current optimal X and
C                         the step length VM. Note that there are
C                         N*NS*NT function evalutations before each
C                         temperature reduction. Finally, notice is
C                         is also given upon achieveing the termination
C                         criteria.
C                     2 - Each new step length (VM), the current optimal
C                         X (XOPT) and the current trial X (X). This
C                         gives the user some idea about how far X
C                         strays from XOPT as well as how VM is adapting
C                         to the function.
C                     3 - Each function evaluation, its acceptance or
C                         rejection and new optima. For many problems,
C                         this option will likely require a small tree
C                         if hard copy is used. This option is best
C                         used to learn about the algorithm. A small
C                         value for MAXEVL is thus recommended when
C                         using IPRINT = 3.
C             Suggested value: 1
C             Note: For a given value of IPRINT, the lower valued
C                   options (other than 0) are utilized.
C    ISEED - The first seed for the random number generator RANMAR.
C             0 .LE. ISEED .LE. 31328. (INT)
C    ISEED2 - The second seed for the random number generator RANMAR.
C             0 .LE. ISEED2 .LE. 30081. Different values for ISEED
C             and ISEED2 will lead to an entirely different sequence
C             of trial points and decisions on downhill moves (when
C             maximizing). See Goffe et al. on how this can be used
C             to test the results of SA. (INT)
C
C  Input/Output Parameters:
C    T - On input, the initial temperature. See Goffe et al. for advice.
C        On output, the final temperature. (DP)
C    VM - The step length vector. On input it should encompass the
C         region of interest given the starting value X. For point
C         X(I), the next trial point is selected is from X(I) - VM(I)
C         to  X(I) + VM(I). Since VM is adjusted so that about half
C         of all points are accepted, the input value is not very
C         important (i.e. if the value is off, SA adjusts VM to the
C         correct value). (DP(N))
C
C  Output Parameters:
C    XOPT - The variables that optimize the function. (DP(N))
C    FOPT - The optimal value of the function. (DP)
C    NACC - The number of accepted function evaluations. (INT)
C    NFCNEV - The total number of function evaluations. In a minor
C             point, note that the first evaluation is not used in the
C             core of the algorithm; it simply initializes the
C             algorithm. (INT).
C    NOBDS - The total number of trial function evaluations that
C            would have been out of bounds of LB and UB. Note that
C            a trial point is randomly selected between LB and UB.
C            (INT)
C    IER - The error return number. (INT)
C          Values: 0 - Normal return; termination criteria achieved.
C                  1 - Number of function evaluations (NFCNEV) is
C                      greater than the maximum number (MAXEVL).
C                  2 - The starting value (X) is not inside the
C                      bounds (LB and UB).
C                  3 - The initial temperature is not positive.
C                  99 - Should not be seen; only used internally.
C
C  Work arrays that must be dimensioned in the calling routine:
C       RWK1 (DP(NEPS))  (FSTAR in SA)
C       RWK2 (DP(N))     (XP    "  " )
C       IWK  (INT(N))    (NACP  "  " )
C
C  Required Functions (included):
C    EXPREP - Replaces the function EXP to avoid under- and overflows.
C             It may have to be modified for non IBM-type main-
C             frames. (DP)
C    RMARIN - Initializes the random number generator RANMAR.
C    RANMAR - The actual random number generator. Note that
C             RMARIN must run first (SA does this). It produces uniform
C             random numbers on [0,1]. These routines are from
C             Usenet's comp.lang.fortran. For a reference, see
C             "Toward a Universal Random Number Generator"
C             by George Marsaglia and Arif Zaman, Florida State
C             University Report: FSU-SCRI-87-50 (1987).
C             It was later modified by F. James and published in
C             "A Review of Pseudo-random Number Generators." For
C             further information, contact stuart@ads.com. These
C             routines are designed to be portable on any machine
C             with a 24-bit or more mantissa. I have found it produces
C             identical results on a IBM 3081 and a Cray Y-MP.
C
C  Required Subroutines (included):
C    PRTVEC - Prints vectors.
C    PRT1 ... PRT10 - Prints intermediate output.
C    FCN - Function to be optimized. The form is
C            SUBROUTINE FCN(N,X,F)
C            INTEGER N
C            DOUBLE PRECISION  X(N), F
C            ...
C            function code with F = F(X)
C            ...
C            RETURN
C            END
C          Note: This is the same form used in the multivariable
C          minimization algorithms in the IMSL edition 10 library.
C
C  Machine Specific Features:
C    1. EXPREP may have to be modified if used on non-IBM type main-
C       frames. Watch for under- and overflows in EXPREP.
C    2. Some FORMAT statements use G25.18; this may be excessive for
C       some machines.
C    3. RMARIN and RANMAR are designed to be protable; they should not
C       cause any problems.

      SUBROUTINE SA(N,X,MAX,RT,EPS,NS,NT,NEPS,MAXEVL,LB,UB,C,IPRINT,
     1              ISEED,ISEED2,T,VM,XOPT,FOPT,NACC,NFCNEV,NOBDS,IER,
     2              FSTAR,XP,NACP)

C  Type all external variables.
      DOUBLE PRECISION  X(N), LB(N), UB(N), C(N), VM(N), FSTAR(NEPS),
     1                  XOPT(N), XP(N), T, EPS, RT, FOPT
      INTEGER  NACP(N), N, NS, NT, NEPS, NACC, MAXEVL, IPRINT,
     1         NOBDS, IER, NFCNEV, ISEED, ISEED2
      LOGICAL  MAX

C  Type all internal variables.
      DOUBLE PRECISION  F, FP, P, PP, RATIO
      INTEGER  NUP, NDOWN, NREJ, NNEW, LNOBDS, H, I, J, M
      LOGICAL  QUIT

C  Type all functions.
      DOUBLE PRECISION  EXPREP
      REAL  RANMAR

C  Initialize the random number generator RANMAR.
      CALL RMARIN(ISEED,ISEED2)

C  Set initial values.
      NACC = 0
      NOBDS = 0
      NFCNEV = 0
      IER = 99

      DO 10, I = 1, N
         XOPT(I) = X(I)
         NACP(I) = 0
10    CONTINUE

      DO 20, I = 1, NEPS
         FSTAR(I) = 1.0D+20
20    CONTINUE 

C  If the initial temperature is not positive, notify the user and 
C  return to the calling routine. 
      IF (T .LE. 0.0) THEN
         WRITE(*,'(/,''  THE INITIAL TEMPERATURE IS NOT POSITIVE. ''
     1             /,''  RESET THE VARIABLE T. ''/)')
         IER = 3
         RETURN
      END IF

C  If the initial value is out of bounds, notify the user and return
C  to the calling routine.
      DO 30, I = 1, N
         IF ((X(I) .GT. UB(I)) .OR. (X(I) .LT. LB(I))) THEN
            CALL PRT1
            IER = 2
            RETURN
         END IF
30    CONTINUE

C  Evaluate the function with input X and return value as F.
      CALL FCN(N,X,F)

C  If the function is to be minimized, switch the sign of the function.
C  Note that all intermediate and final output switches the sign back
C  to eliminate any possible confusion for the user.
      IF(.NOT. MAX) F = -F
      NFCNEV = NFCNEV + 1
      FOPT = F
      FSTAR(1) = F
      IF(IPRINT .GE. 1) CALL PRT2(MAX,N,X,F)

C  Start the main loop. Note that it terminates if (i) the algorithm
C  succesfully optimizes the function or (ii) there are too many
C  function evaluations (more than MAXEVL).
100   NUP = 0
      NREJ = 0
      NNEW = 0
      NDOWN = 0
      LNOBDS = 0

      DO 400, M = 1, NT
         DO 300, J = 1, NS
            DO 200, H = 1, N

C  Generate XP, the trial value of X. Note use of VM to choose XP.
               DO 110, I = 1, N
                  IF (I .EQ. H) THEN
                     XP(I) = X(I) + (RANMAR()*2.- 1.) * VM(I)
                  ELSE
                     XP(I) = X(I)
                  END IF

C  If XP is out of bounds, select a point in bounds for the trial.
                  IF((XP(I) .LT. LB(I)) .OR. (XP(I) .GT. UB(I))) THEN
                    XP(I) = LB(I) + (UB(I) - LB(I))*RANMAR()
                    LNOBDS = LNOBDS + 1
                    NOBDS = NOBDS + 1
                    IF(IPRINT .GE. 3) CALL PRT3(MAX,N,XP,X,FP,F)
                  END IF
110            CONTINUE

C  Evaluate the function with the trial point XP and return as FP.
               CALL FCN(N,XP,FP)

               IF(.NOT. MAX) FP = -FP
               NFCNEV = NFCNEV + 1
               IF(IPRINT .GE. 3) CALL PRT4(MAX,N,XP,X,FP,F)

C  If too many function evaluations occur, terminate the algorithm.
               IF(NFCNEV .GE. MAXEVL) THEN
                  CALL PRT5
                  IF (.NOT. MAX) FOPT = -FOPT
                  IER = 1
                  RETURN
               END IF

C  Accept the new point if the function value increases.
               IF(FP .GE. F) THEN
                  IF(IPRINT .GE. 3) THEN
                     WRITE(*,'(''  POINT ACCEPTED'')')
                  END IF
                  DO 120, I = 1, N
                     X(I) = XP(I)
120               CONTINUE
                  F = FP
                  NACC = NACC + 1
                  NACP(H) = NACP(H) + 1
                  NUP = NUP + 1

C  If greater than any other point, record as new optimum.
                  IF (FP .GT. FOPT) THEN
                  IF(IPRINT .GE. 3) THEN
                     WRITE(*,'(''  NEW OPTIMUM'')')
                  END IF
                     DO 130, I = 1, N
                        XOPT(I) = XP(I)
130                  CONTINUE
                     FOPT = FP
                     NNEW = NNEW + 1
                  END IF

C  If the point is lower, use the Metropolis criteria to decide on
C  acceptance or rejection.
               ELSE
                  P = EXPREP((FP - F)/T)
                  PP = RANMAR()
                  IF (PP .LT. P) THEN
                     IF(IPRINT .GE. 3) CALL PRT6(MAX)
                     DO 140, I = 1, N
                        X(I) = XP(I)
140                  CONTINUE
                     F = FP
                     NACC = NACC + 1
                     NACP(H) = NACP(H) + 1
                     NDOWN = NDOWN + 1
                  ELSE
                     NREJ = NREJ + 1
                     IF(IPRINT .GE. 3) CALL PRT7(MAX)
                  END IF
               END IF

200         CONTINUE
300      CONTINUE

C  Adjust VM so that approximately half of all evaluations are accepted.
         DO 310, I = 1, N
            RATIO = DFLOAT(NACP(I)) /DFLOAT(NS)
            IF (RATIO .GT. .6) THEN
               VM(I) = VM(I)*(1. + C(I)*(RATIO - .6)/.4)
            ELSE IF (RATIO .LT. .4) THEN
               VM(I) = VM(I)/(1. + C(I)*((.4 - RATIO)/.4))
            END IF
            IF (VM(I) .GT. 0.5d0*(UB(I)-LB(I))) THEN
               VM(I) = 0.5d0*(UB(I) - LB(I))
            END IF
310      CONTINUE

         IF(IPRINT .GE. 2) THEN
            CALL PRT8(N,VM,XOPT,X)
         END IF

         DO 320, I = 1, N
            NACP(I) = 0
320      CONTINUE

400   CONTINUE

      IF(IPRINT .GE. 1) THEN
         CALL PRT9(MAX,N,T,XOPT,VM,FOPT,NUP,NDOWN,NREJ,LNOBDS,NNEW)
      END IF

C  Check termination criteria.
      QUIT = .FALSE.
      FSTAR(1) = F
      IF ((FOPT - FSTAR(1)) .LE. EPS) QUIT = .TRUE.
      DO 410, I = 1, NEPS
         IF (ABS(F - FSTAR(I)) .GT. EPS) QUIT = .FALSE.
410   CONTINUE

C  Terminate SA if appropriate.
      IF (QUIT) THEN
         DO 420, I = 1, N
            X(I) = XOPT(I)
420      CONTINUE
         IER = 0
         IF (.NOT. MAX) FOPT = -FOPT
         IF(IPRINT .GE. 1) CALL PRT10
         RETURN
      END IF

C  If termination criteria is not met, prepare for another loop.
      NT = NT-5
      if (nt.eq.0) then
         write(*,'("***** NT reached 0 ******")')
         if (.not. max) fopt = -fopt
         return
      endif
      T = RT*T
      DO 430, I = NEPS, 2, -1
         FSTAR(I) = FSTAR(I-1)
430   CONTINUE
      F = FOPT
      DO 440, I = 1, N
         X(I) = XOPT(I)
440   CONTINUE

C  Loop again.
      GO TO 100

      END

      FUNCTION  EXPREP(RDUM)
C  This function replaces exp to avoid under- and overflows and is
C  designed for IBM 370 type machines. It may be necessary to modify
C  it for other machines. Note that the maximum and minimum values of
C  EXPREP are such that they has no effect on the algorithm.

      DOUBLE PRECISION  RDUM, EXPREP

      IF (RDUM .GT. 174.) THEN
         EXPREP = 3.69D+75
      ELSE IF (RDUM .LT. -180.) THEN
         EXPREP = 0.0
      ELSE
         EXPREP = EXP(RDUM)
      END IF

      RETURN
      END

      subroutine RMARIN(IJ,KL)
C  This subroutine and the next function generate random numbers. See
C  the comments for SA for more information. The only changes from the
C  orginal code is that (1) the test to make sure that RMARIN runs first
C  was taken out since SA assures that this is done (this test didn't
C  compile under IBM's VS Fortran) and (2) typing ivec as integer was
C  taken out since ivec isn't used. With these exceptions, all following
C  lines are original.

C This is the initialization routine for the random number generator
C     RANMAR()
C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
C                                                      0 <= KL <= 30081
      real U(97), C, CD, CM
      integer I97, J97
      common /raset1/ U, C, CD, CM, I97, J97
      if( IJ .lt. 0  .or.  IJ .gt. 31328  .or.
     *    KL .lt. 0  .or.  KL .gt. 30081 ) then
          print '(A)', ' The first random number seed must have a value
     *between 0 and 31328'
          print '(A)',' The second seed must have a value between 0 and
     *30081'
            stop
      endif
      i = mod(IJ/177, 177) + 2
      j = mod(IJ    , 177) + 2
      k = mod(KL/169, 178) + 1
      l = mod(KL,     169)
      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
3        continue
         U(ii) = s
2     continue
      C = 362436.0 / 16777216.0
      CD = 7654321.0 / 16777216.0
      CM = 16777213.0 /16777216.0
      I97 = 97
      J97 = 33
      return
      end

      function ranmar()
      real U(97), C, CD, CM
      integer I97, J97
      common /raset1/ U, C, CD, CM, I97, J97
         uni = U(I97) - U(J97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         U(I97) = uni
         I97 = I97 - 1
         if(I97 .eq. 0) I97 = 97
         J97 = J97 - 1
         if(J97 .eq. 0) J97 = 97
         C = C - CD
         if( C .lt. 0.0 ) C = C + CM
         uni = uni - C
         if( uni .lt. 0.0 ) uni = uni + 1.0
         RANMAR = uni
      return
      END

      SUBROUTINE PRT1
C  This subroutine prints intermediate output, as does PRT2 through
C  PRT10. Note that if SA is minimizing the function, the sign of the
C  function value and the directions (up/down) are reversed in all
C  output to correspond with the actual function optimization. This
C  correction is because SA was written to maximize functions and
C  it minimizes by maximizing the negative a function.

      WRITE(*,'(/,''  THE STARTING VALUE (X) IS OUTSIDE THE BOUNDS ''
     1          /,''  (LB AND UB). EXECUTION TERMINATED WITHOUT ANY''
     2          /,''  OPTIMIZATION. RESPECIFY X, UB OR LB SO THAT  ''
     3          /,''  LB(I) .LT. X(I) .LT. UB(I), I = 1, N. ''/)')

      RETURN
      END

      SUBROUTINE PRT2(MAX,N,X,F)

      DOUBLE PRECISION  X(*), F
      INTEGER  N
      LOGICAL  MAX

c      WRITE(*,'(''  '')')
      CALL PRTVEC(X,N,'INITIAL X')
      IF (MAX) THEN
         WRITE(*,'(/''  INITIAL F:  '', G25.18)') F
      ELSE
         WRITE(*,'(/''  INITIAL F:  '', G25.18)') -F
      END IF

      RETURN
      END

      SUBROUTINE PRT3(MAX,N,XP,X,FP,F)

      DOUBLE PRECISION  XP(*), X(*), FP, F
      INTEGER  N
      LOGICAL  MAX

      WRITE(*,'(''  '')')
      CALL PRTVEC(X,N,'CURRENT X')
      IF (MAX) THEN
         WRITE(*,'(''  CURRENT F: '',G25.18)') F
      ELSE
         WRITE(*,'(''  CURRENT F: '',G25.18)') -F
      END IF
      CALL PRTVEC(XP,N,'TRIAL X')
      WRITE(*,'(''  POINT REJECTED SINCE OUT OF BOUNDS'')')

      RETURN
      END

      SUBROUTINE PRT4(MAX,N,XP,X,FP,F)

      DOUBLE PRECISION  XP(*), X(*), FP, F
      INTEGER  N
      LOGICAL  MAX

      WRITE(*,'(''  '')')
      CALL PRTVEC(X,N,'CURRENT X')
      IF (MAX) THEN
         WRITE(*,'(''  CURRENT F: '',G25.18)') F
         CALL PRTVEC(XP,N,'TRIAL X')
         WRITE(*,'(''  RESULTING F: '',G25.18)') FP
      ELSE
         WRITE(*,'(''  CURRENT F: '',G25.18)') -F
         CALL PRTVEC(XP,N,'TRIAL X')
         WRITE(*,'(''  RESULTING F: '',G25.18)') -FP
      END IF

      RETURN
      END

      SUBROUTINE PRT5

      WRITE(*,'(/,''  TOO MANY FUNCTION EVALUATIONS; CONSIDER ''
     1          /,''  INCREASING MAXEVL OR EPS, OR DECREASING ''
     2          /,''  NT OR RT. THESE RESULTS ARE LIKELY TO BE ''
     3          /,''  POOR.'',/)')

      RETURN
      END

      SUBROUTINE PRT6(MAX)

      LOGICAL  MAX

      IF (MAX) THEN
         WRITE(*,'(''  THOUGH LOWER, POINT ACCEPTED'')')
      ELSE
         WRITE(*,'(''  THOUGH HIGHER, POINT ACCEPTED'')')
      END IF

      RETURN
      END

      SUBROUTINE PRT7(MAX)

      LOGICAL  MAX

      IF (MAX) THEN
         WRITE(*,'(''  LOWER POINT REJECTED'')')
      ELSE
         WRITE(*,'(''  HIGHER POINT REJECTED'')')
      END IF

      RETURN
      END

      SUBROUTINE PRT8(N,VM,XOPT,X)

      DOUBLE PRECISION  VM(*), XOPT(*), X(*)
      INTEGER  N

      WRITE(*,'(/,
     1  '' INTERMEDIATE RESULTS AFTER STEP LENGTH ADJUSTMENT'',/)')
      CALL PRTVEC(VM,N,'NEW STEP LENGTH (VM)')
      CALL PRTVEC(XOPT,N,'CURRENT OPTIMAL X')
      CALL PRTVEC(X,N,'CURRENT X')
      WRITE(*,'('' '')')

      RETURN
      END

      SUBROUTINE PRT9(MAX,N,T,XOPT,VM,FOPT,NUP,NDOWN,NREJ,LNOBDS,NNEW)

      DOUBLE PRECISION  XOPT(*), VM(*), T, FOPT
      INTEGER  N, NUP, NDOWN, NREJ, LNOBDS, NNEW, TOTMOV
      LOGICAL  MAX

      TOTMOV = NUP + NDOWN + NREJ

      WRITE(*,'(/,
     1  '' INTERMEDIATE RESULTS BEFORE NEXT TEMPERATURE REDUCTION'',/)')
      WRITE(*,'(''  CURRENT TEMPERATURE:            '',G12.5)') T
      IF (MAX) THEN
         WRITE(*,'(''  MAX FUNCTION VALUE SO FAR:  '',G25.18)') FOPT
         WRITE(*,'(''  TOTAL MOVES:                '',I8)') TOTMOV
         WRITE(*,'(''     UPHILL:                  '',I8)') NUP
         WRITE(*,'(''     ACCEPTED DOWNHILL:       '',I8)') NDOWN
         WRITE(*,'(''     REJECTED DOWNHILL:       '',I8)') NREJ
         WRITE(*,'(''  OUT OF BOUNDS TRIALS:       '',I8)') LNOBDS
         WRITE(*,'(''  NEW MAXIMA THIS TEMPERATURE:'',I8)') NNEW
      ELSE
         WRITE(*,'(''  MIN FUNCTION VALUE SO FAR:  '',G25.18)') -FOPT
         WRITE(*,'(''  TOTAL MOVES:                '',I8)') TOTMOV
         WRITE(*,'(''     DOWNHILL:                '',I8)')  NUP
         WRITE(*,'(''     ACCEPTED UPHILL:         '',I8)')  NDOWN
         WRITE(*,'(''     REJECTED UPHILL:         '',I8)')  NREJ
         WRITE(*,'(''  TRIALS OUT OF BOUNDS:       '',I8)')  LNOBDS
         WRITE(*,'(''  NEW MINIMA THIS TEMPERATURE:'',I8)')  NNEW
      END IF
      CALL PRTVEC(XOPT,N,'CURRENT OPTIMAL X')
      CALL PRTVEC(VM,N,'STEP LENGTH (VM)')
      WRITE(*,'('' '')')

      RETURN
      END

      SUBROUTINE PRT10

      WRITE(*,'(/,''  SA ACHIEVED TERMINATION CRITERIA. IER = 0. '',
     &/)')

      RETURN
      END

      SUBROUTINE PRTVEC(VECTOR,NCOLS,NAME)
C  This subroutine prints the double precision vector named VECTOR.
C  Elements 1 thru NCOLS will be printed. NAME is a character variable
C  that describes VECTOR. Note that if NAME is given in the call to
C  PRTVEC, it must be enclosed in quotes. If there are more than 10
C  elements in VECTOR, 10 elements will be printed on each line.

      INTEGER NCOLS
      DOUBLE PRECISION VECTOR(NCOLS)
      CHARACTER *(*) NAME

      WRITE(*,1001) NAME

      IF (NCOLS .GT. 6) THEN
         LINES = INT(NCOLS/6.)

         DO 100, I = 1, LINES
            LL = 6*(I - 1)
            WRITE(*,1000) (VECTOR(J),J = 1+LL, 6+LL)
  100    CONTINUE

         WRITE(*,1000) (VECTOR(J),J = 7+LL, NCOLS)
      ELSE
         WRITE(*,1000) (VECTOR(J),J = 1, NCOLS)
      END IF

 1000 FORMAT( 10(G12.5,1X))
 1001 FORMAT(/,25X,A)

      RETURN
      END

